#include <iostream>
#include <complex>
#include "nnec/nnec_kernels.h"

using namespace std;
using namespace NNEC;

int main()
{
    int type;
    double rho, z, d, a, z1, z2, step;
    double ea, eA, eb, eB;
    double v_abs;

    complex<double> g, ga, gA, gb, gB;

    complex<double> (* pG)(double rho, double a, double z, double z0);
    complex<double> (* pGA)(double rho, double a, double z, double z0);
    complex<double> (* pGa)(double rho, double z, double z0);
    complex<double> (* pGB)(double rho, double a, double z, double z0);
    complex<double> (* pGb)(double rho, double a, double z, double z0);

    d = 0.005;
    a = 0.001;
    type = 1;

    switch(type) {
        case 1:
            pG = G1;
            pGA = G1A;
            pGa = G1a;
            pGB = G1B;
            pGb = G1b;
            break;
        case 2:
            pG = G2;
            pGA = G2A;
            pGa = G2a;
            pGB = G2B;
            pGb = G2b;
            break;
        case 3:
            pG = G3;
            pGA = G3A;
            pGa = G3a;
            pGB = G3B;
            pGb = G3b;
            break;
        default:
            break;
    }

    if(type >= 1 && type <= 3) {
        rho = 0.5 * a;
        z1 = 0.0;
        z2 = 0.0;
        step = 0.01 * a;
        while(step <= 20e-2) {
            /*if(abs(step - a) / a < 1e-3 ) {
                rho = (1 + 1e-3) * a;
            }
            else {
                rho = step;
            }*/
            z = step;

            g = pG(rho, a, z, z1);
            gA = pGA(rho, a, z, z1);
            ga = pGa(rho, z, z1);
            gB = pGB(rho, a, z, z1);
            gb = pGb(rho, a, z, z1);

            v_abs = abs(g);
            if(v_abs != 0.0) {
                eA = abs(gA - g) / v_abs;
                ea = abs(ga - g) / v_abs;
                eB = abs(gB - g) / v_abs;
                eb = abs(gb - g) / v_abs;
            }
            else {
                eA = 0.0;
                ea = 0.0;
                eB = 0.0;
                eb = 0.0;
            }

            cout.precision(16);
            //cout << rho << " ";
            cout << z << " ";
            cout << real(g) << " " << imag(g) << " ";
            cout << real(gA) << " " << imag(gA) << " ";
            cout << real(ga) << " " << imag(ga) << " ";
            cout << real(gB) << " " << imag(gB) << " ";
            cout << real(gb) << " " << imag(gb) << " ";
            cout << eA << " " << ea << " " << eB << " " << eb << endl;

            step = step + 0.01 * a;
        }
    }
    else if(type == 4) {
        z1 = -d;
        z2 = d;
        rho = 3.0 * a;
        step = 0.01 * a;
        while(step <= 20e-3) {
            /*if(abs(step - a) / a < 1e-3 ) {
                rho = (1 + 1e-3) * a;
            }
            else {
                rho = step;
            }*/
            z = step;
            g = G4(rho, a, z, z1, z2);
            gA = G4A(rho, a, z, z1, z2);
            ga = G4a(rho, z, z1, z2);
            gB = complex<double>(0.0, 0.0);
            gb = G4b(rho, a, z, z1, z2);

            v_abs = abs(g);
            if(v_abs != 0.0) {
                eA = abs(gA - g) / v_abs;
                ea = abs(ga - g) / v_abs;
                eB = abs(gB - g) / v_abs;
                eb = abs(gb - g) / v_abs;
            }
            else {
                eA = 0.0;
                ea = 0.0;
                eB = 0.0;
                eb = 0.0;
            }

            cout.precision(16);
            //cout << rho << " ";
            cout << z << " ";
            cout << real(g) << " " << imag(g) << " ";
            cout << real(gA) << " " << imag(gA) << " ";
            cout << real(ga) << " " << imag(ga) << " ";
            cout << real(gB) << " " << imag(gB) << " ";
            cout << real(gb) << " " << imag(gb) << " ";
            cout << eA << " " << ea << " " << eB << " " << eb << endl;

            step = step + 0.01 * a;
        }
    }

    return 0;
}